#!/usr/bin/env python3

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import cartopy.crs as ccrs
import sys

lon = 0
lat = 0
proj = ccrs.Orthographic(central_longitude=lon, central_latitude=lat)
# proj = ccrs.PlateCarree()
proj._threshold /= 100

plt.figure(figsize=(8, 8))
ax = plt.axes(projection=proj)
ax.coastlines(resolution='110m', lw=0.1)
ax.set_global()

mesh = xr.open_dataset(sys.argv[1])

ax.scatter(mesh.lon, mesh.lat, s=0.5, c='r', transform=ccrs.PlateCarree())
ax.scatter(mesh.lon_eq, mesh.lat_eq, s=0.5, c='g', transform=ccrs.PlateCarree())
ax.scatter(mesh.lon_cq, mesh.lat_cq, s=0.5, c='g', transform=ccrs.PlateCarree())

norm = mpl.colors.Normalize(vmin=0.5, vmax=1)
cmap = cm.rainbow
color = cm.ScalarMappable(norm=norm, cmap=cmap)
for k in range(mesh.dom.size):
	for j in range(mesh.y.size):
		for i in range(mesh.x.size):
			lon = list(mesh.lon_v[k,j,i,:]); lon.append(mesh.lon_v[k,j,i,0])
			lat = list(mesh.lat_v[k,j,i,:]); lat.append(mesh.lat_v[k,j,i,0])
			ax.plot(lon, lat, c='b', lw=0.3, transform=ccrs.Geodetic())
			# ax.fill(lon, lat, c=color.to_rgba(mesh.J_c[k,j,i]), transform=ccrs.Geodetic())

# Plot domain outlines.
#rot_proj = ccrs.Orthographic(central_longitude=mesh.rlon, central_latitude=mesh.rlat) # Not used yet!
outlines = np.array([
	[[315, -35.2643896827547], [45 , -35.2643896827547]],
	[[45 , -35.2643896827547], [45 ,  35.2643896827547]],
	[[45 ,  35.2643896827547], [315,  35.2643896827547]],
	[[315,  35.2643896827547], [315, -35.2643896827547]],
	[[45 , -35.2643896827547], [135, -35.2643896827547]],
	[[135, -35.2643896827547], [135,  35.2643896827547]],
	[[135,  35.2643896827547], [45 ,  35.2643896827547]],
	[[135, -35.2643896827547], [225, -35.2643896827547]],
	[[225, -35.2643896827547], [225,  35.2643896827547]],
	[[225,  35.2643896827547], [135,  35.2643896827547]],
	[[225, -35.2643896827547], [315, -35.2643896827547]],
	[[315,  35.2643896827547], [225,  35.2643896827547]]
])
for line in outlines:
	ax.plot(line[:,0], line[:,1], c='black', lw=0.6, transform=ccrs.Geodetic())

plt.show()
